Determining if particles are spirally arranged

The goal was to determine if particles on the electron microscope images follow spiral or concetric arrangement.

import matplotlib as mpl
from matplotlib import pyplot as plt
import numpy as np
import cv2
import os
import copy
from scipy.optimize import curve_fit

import nanoparticles as nano
import importlib
importlib.reload(nano)

<module 'nanoparticles' from 'h:\\My Drive\\Projekty\\Small projects\\spiral\\nanoparticles.py'>

Prepare data

Load and preprocess images

file = 'W14_J1_3_47 [1].tif'
img = cv2.imread('data/' + file)
mask = cv2.imread('data/' + file.replace('.', '_mask.'))
mask = cv2.cvtColor(mask, cv2.COLOR_RGB2GRAY)
centers = cv2.imread('data/' + file.replace('.', '_centers.'))

img = cv2.flip(img, 1)
mask = cv2.flip(mask, 1)
centers = cv2.flip(centers, 1)
imgGray = cv2.cvtColor(img.copy(), cv2.COLOR_RGB2GRAY)
scale, img_scale = nano.find_scalebar(
    imgGray.copy(), 200, showscale=True)  # scale is in nm/px


print(f'Scale is {scale:.4} nm/px')
nano.plot_images([img_scale, mask, centers], names=['Raw Image (+scalebar)',
                 'Manually selected regions', 'Manually selected centers of spirals'], size=14)

Scale is 0.3035 nm/px

def get_spirals_centers(centers):
    centersGray = cv2.cvtColor(centers, cv2.COLOR_RGB2GRAY)
    ret, thrCen = cv2.threshold(255-centersGray, 5, 255, cv2.THRESH_BINARY)
    contours, hierarchy = cv2.findContours(
        thrCen, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

    spirals = []
    # colored = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
    for c in contours:
        # compute the center of the contour
        M = cv2.moments(c)
        cX = int(M["m10"] / M["m00"])
        cY = int(M["m01"] / M["m00"])
        spirals.append([cX, cY])

    return spirals


spirals = get_spirals_centers(centers)
fewSpirals = [nano.get_aoi(img, spiral) for spiral in spirals[:5]]
nano.plot_images(fewSpirals, size=8, rowsplit=5)

Thresholding and finding particles

blur = cv2.medianBlur(imgGray, 9)

# use adaptive gaussian thresholding
blocksize = 251
thr = cv2.adaptiveThreshold(
    blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, blocksize, 2)

thr = np.where(thr == 255, 0, 255).astype(np.uint8)

# erode particles
kernel = np.ones((5, 5), np.uint8)
thrEroded = cv2.erode(thr.copy(), kernel, iterations=1)

# select particles
dMin, dMax = 1, 7  # min and max diameter of particles in nm
gCon, bCon = nano.find_contours(
    thrEroded, scale, areaMin=np.pi*(dMin/2)**2, areaMax=np.pi*(dMax/2)**2, circMin=0.5)
rgbColors = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
rgbColors = cv2.drawContours(rgbColors, gCon, -1, (0,255,0), 1)
rgbColors = cv2.drawContours(rgbColors, bCon, -1, (255,0,0), 1)
particles = nano.get_centers(gCon)

names = ['Blurred', 'Threshold', 'Eroded', 'Particles']
nano.plot_images([nano.get_aoi(image, spirals[1], size = 500) for image in [blur, thr, thrEroded, rgbColors]], names = names, size = 12)

Analyze particles

# filter particles that are outside mask
spiralPrtcs = []
for p in particles:
    if mask[p[1], p[0]] != 0:
        spiralPrtcs.append(p)        
spiralPrtcs = np.array(spiralPrtcs)

# assign to spiral based on distance to centers
particlesAssigned = [[] for _ in range(len(spirals))]
for p in spiralPrtcs:
    distances = []
    for s in spirals:
        distances.append(nano.point_distance(p, s))
        
    minD = np.argmin(np.array(distances))
    particlesAssigned[minD].append(p)

plt.imshow(imgGray, cmap = 'gray')

for prt in particlesAssigned:
    plt.scatter(np.array(prt)[:, 0], np.array(prt)[:, 1], s = 0.3)
    
plt.scatter(np.array(spirals)[:, 0], np.array(spirals)[:, 1], c = 'black', s = 15)
for i, spiral in enumerate(spirals):
    x, y = spiral[0], spiral[1]
    plt.text(x, y, str(i), color = 'black')

rgbColors = cv2.cvtColor(imgGray,cv2.COLOR_GRAY2RGB)
for n in range(len(spirals)):
    color = nano.random_color()
    for particle in particlesAssigned[n]:
        rgbColors = cv2.circle(rgbColors, particle, 5, color, -1)
nano.plot_images([nano.get_aoi(rgbColors, spiral) for spiral in spirals[:4]], size = 12)

# add polar coordinates
particlesDetails = []
for i in range(len(spirals)):
    spiral = spirals[i]
    particles = []
    for x, y in particlesAssigned[i]:
        dist, alfa = nano.cart2pol(x-spiral[0], y-spiral[1])
        particles.append([x, y, dist, alfa])
    particlesDetails.append(np.array(particles))

Assign particles to rows

We are looking for rows of particles that follow along the curve around the center of spiral. We first select one particle (that is closes to center) and then find next particle like this:

How alfa was calculated

maxDist = 3.5/scale
idealSpiralDist = 9/scale
drs = []

# grouping particles into rows that follow curve around the center of spiral
def group_particles(particlesDetails):
    particlesGrouped = []
    for i in range(len(particlesDetails)):
        group = []
        particles = copy.copy(particlesDetails[i])
        particles = np.array(particles)
        particles = list(particles[particles[:, 2].argsort()])  # sort by distance to center
        
        prt = particles[0]  # get first particle
        row = [prt]  # creates initial row
        del particles[0]
        direction = 0  # direction (clockwise or anti) should be unifrom within row of particles
        while len(particles) > 0:
            x, y, r, theta = row[-1]
            # alfa is the expect angle of ideal neighbouring particle at that particular distance to center (r)
            alfa = 2*np.arcsin(idealSpiralDist/(2*r))  
            
            # if direction is selected
            if direction != 0:
                x2, y2 = nano.pol2cart(r, theta+alfa*direction)
                dsts = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2, y2)) for p in particles]
                minD = np.argmin(dsts)
                if dsts[minD] < maxDist:
                    row.append(particles[minD])
                    del particles[minD]
                else:
                    if len(row) > 1: group.append(row)  # only take rows with more than 1 particle
                    row = [particles[0]]  # starts new row
                    del particles[0]
                    direction = 0
                    
            else:  # find direction (check on which side best particle match is closest to ideal)
                # positive direction
                x2P, y2P = nano.pol2cart(r, theta+alfa)
                dstsP = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2P, y2P)) for p in particles]
                minDP = np.argmin(dstsP)

                # negative direction
                x2N, y2N = nano.pol2cart(r, theta-alfa)
                dstsN = [nano.point_distance((p[0] - spirals[i][0], p[1] - spirals[i][1]), (x2N, y2N)) for p in particles]
                minDN = np.argmin(dstsN)

                # select
                if min(dstsN[minDN], dstsP[minDP]) < maxDist:
                    if dstsN[minDN] < dstsP[minDP]:
                        row.append(particles[minDN])
                        del particles[minDN]
                        direction = -1
                    else:
                        row.append(particles[minDP])
                        del particles[minDP]
                        direction = 1
                    drs.append(direction)
                else:
                    row = [particles[0]]
                    del particles[0]
                    direction = 0
                    
        particlesGrouped.append(group)
    return particlesGrouped


particlesGrouped = group_particles(particlesDetails)

C:\Users\szust\AppData\Local\Temp\ipykernel_23676\666208268.py:21: RuntimeWarning: invalid value encountered in arcsin
  alfa = 2*np.arcsin(idealSpiralDist/(2*r))

# just drawing
work = copy.copy(img)

for n, spiral in enumerate(particlesGrouped):
    work = cv2.circle(work, spirals[n], 10, (255, 255, 255), -1)
    for group in spiral:
        color = nano.random_color()
        for i in range(0, len(group)-1):
            p1 = np.array(group[i][:2])
            p2 = np.array(group[i+1][:2])
            work = cv2.line(work, p1.astype(int), p2.astype(int), color, 2)
            work = cv2.circle(work, p1.astype(int), 4, color, -1)
        work = cv2.circle(work, p2.astype(int), 4, color, -1)
            
nano.plot_images([nano.get_aoi(work, spiral) for spiral in spirals[:4]], size = 12)

Analyzing particle arrangement

def linear_func(x, a, b):
    return a*x + b


def fix_angles_inrow(angles):
    '''
    add/substract 2*np.pi from list of angles after going through 0/2*np.pi (idicated when theta difference is > np.pi)
    '''
    newAngles = [angles[0]]
    adder = 0
    for i in range(1, len(angles)):
        if angles[i-1] - angles[i] > np.pi:
            adder += 2*np.pi
        if angles[i] - angles[i-1] > np.pi:
            adder -= 2*np.pi
        newAngles.append(angles[i] + adder)
    return np.array(newAngles)


def calc_and_plot(rows, work, center):
    fig, axs = plt.subplots(ncols = 2, nrows = 1, figsize = (12, 6*1))
    factorRng = 5  # to determine color
    coefs, thetasSelected, distancesSelected = [], [], []
    for row in rows:
        row = np.array(row)
        if len(row) > 2 and min(row[:, 2])*scale > 10:  # take only rows that are not too close to center
            distances = row[:, 2]
            thetas = row[:, 3]
            thetas += np.pi  # to not have negative theta (does not matter really)

            thetas = fix_angles_inrow(thetas)

            # fitting to linear equation to get "a" coeficient (y = ax + b)
            popt, pcov = curve_fit(linear_func, thetas, distances, bounds=(-np.inf, np.inf))
            a = popt[0]*scale
            coefs.append(a)
            meanCoef = np.mean(coefs)

            # draw this row on chart
            c = mpl.cm.get_cmap('coolwarm')((factorRng+a)/(2*factorRng))
            axs[0].text(np.mean(thetas), np.mean(distances)*scale, f'{a*scale:.1f}', fontsize = 7)
            axs[0].scatter(thetas, distances*scale, s = 6, color = c)
            axs[0].plot(thetas, distances*scale, color = c, lw = 1.7)
            axs[0].set_xlabel('Theta [rad]')
            axs[0].set_ylabel('Distance [nm]')
            axs[0].set_title('Fitted a (dist = a * theta + b)')

            # draw this row on image
            ccv = (np.array(c[:3])*255).astype(int)  # color used by cv2
            ccv = tuple(map(int, ccv))
            for ip in range(0, len(row)-1):
                p1 = np.array(row[ip][:2])
                p2 = np.array(row[ip+1][:2])
                work = cv2.line(work, p1.astype(int), p2.astype(int), ccv, 5)
                work = cv2.circle(work, p1.astype(int), int(2/scale), ccv, -2)
            work = cv2.circle(work, p2.astype(int), int(2/scale), ccv, -2)

            # draw colored dot at center of spiral
            c = mpl.cm.get_cmap('coolwarm')((factorRng+meanCoef)/(2*factorRng))
            ccv = (np.array(c[:3])*255).astype(int)  # color used by cv2
            ccv = tuple(map(int, ccv))
            work = cv2.circle(work, center, int(4/scale), (0,0,0), -2)
            work = cv2.circle(work, center, int(3/scale), (ccv), -2)
            
            axs[1].imshow(nano.get_aoi(work, center))
            axs[1].set_title('Corresponding image')
            thetasSelected.append(thetas)
            distancesSelected.append(distances)
    plt.close()
    return fig, [coefs, thetasSelected, distancesSelected]

Result

fig, result = calc_and_plot(particlesGrouped[2], copy.copy(img), spirals[2])
fig

fig, result = calc_and_plot(particlesGrouped[6], copy.copy(img), spirals[6])
fig

fig, result = calc_and_plot(particlesGrouped[8], copy.copy(img), spirals[8])
fig

For all spirals

work = copy.copy(img)
results = []
for i in range(len(spirals)):
    fig, result = calc_and_plot(particlesGrouped[i], work, spirals[i])
    results.append(result)

nano.plot_images([work], size = 30)

fig, axs = plt.subplots(ncols = 2, nrows = 1, figsize = (12, 6*1))
factorRng = 5
allCoefs = []

for result in results:
    coefs, thetasSelected, distancesSelected = copy.copy(result)
    for coef, thetas, distances in zip(coefs, thetasSelected, distancesSelected):
        c = mpl.cm.get_cmap('coolwarm')((factorRng+coef)/(2*factorRng))
        firstPoint = np.argmin(thetas)
        thetas -= thetas[firstPoint]
        distances -= distances[firstPoint]
        axs[0].scatter(thetas, distances*scale, s = 6, color = c)
        axs[0].plot(thetas, distances*scale, color = c, lw = 1.7)
        axs[0].set_xlabel('Theta [rad]')
        axs[0].set_ylabel('ΔR [nm]')
        axs[0].set_title('Fitted coefficient (R = a * theta + b)')
        allCoefs.append(coef)

_ = axs[1].hist(allCoefs, bins = 40)
axs[1].set_title('Coefficient histogram')
axs[1].set_xlabel('Coefficient')
axs[1].set_ylabel('Count')

Text(0, 0.5, 'Count')

CONCLUSION: Depending on the specific spiral, particles can follow right (a > 0, red color) or left (a < 0, blue color) spiral arrangement or be close to concetric (a = 0).